execute mcmc sampling

goMCMC=function(mdl,data,smp=500,wrm=100,th=1){
  mcmc=mdl$sample(
  data=data,
  seed=1,
  chains=4,
  iter_sampling=smp,
  iter_warmup=wrm,
  thin=th,
  refresh=1000
  )
  mcmc
}

see mcmc result and parameters

seeMCMC=function(mcmc,exc='',ch=T){ # exclude 'exc' parameters from seeing
  print(mcmc)
  prs=mcmc$metadata()$model_params[-1] # reject lp__
  for(pr in prs){
    if(grepl('^y',pr)) next # not show predictive value "y*" information
    if(exc!='' && grepl(paste0('^',exc),pr)) next
    drw=mcmc$draws(pr)
    if(ch){
      par(mfrow=c(1,4))
      drw[,1,] |> plot(type='l',main='chain1',ylab=pr)
      drw[,2,] |> plot(type='l',main='chain2',ylab=pr)
      drw[,3,] |> plot(type='l',main='chain3',ylab=pr)
      drw[,4,] |> plot(type='l',main='chain4',ylab=pr)
    }

    par(mfrow=c(1,2))
    drw |> hist(main=pr,xlab='')
    drw |> density() |> plot(main=pr)    
  }
}
fn=function(mdl,data,smp=500,wrm=100){
mle=mdl$optimize(data=data)
print(mle)

mcmc=goMCMC(mdl,data,smp,wrm)

mcmc$metadata()$stan_variables
seeMCMC(mcmc,'m')


y0=mcmc$draws('y0')
smy0=summary(y0)

grid.arrange(
  qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')


m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)

xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)

qplot(x,y,col=I('red'))+
  geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')
}

ex9

explanatory variable obs.x also has noise  
x~N(x0.sx)  
y~N(b0+b1*x0,s)  

ex8-0.stan

data{
  int N;
  vector[N] x;
  vector[N] y;
  int N1;
  vector[N1] x1;
}
parameters{
  real b0;
  real b1;
  real<lower=0> s;
}
model{
  y~normal(b0+b1*x,s);
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0+b1*x[i];
    y0[i]=normal_rng(m0[i],s);
  }
  vector[N1] m1;
  vector[N1] y1;
  for(i in 1:N1){
    m1[i]=b0+b1*x1[i];
    y1[i]=normal_rng(m1[i],s);
  }
}

ex9.stan

data{
  int N;
  vector[N] x;
  vector[N] y;
  int N1;
  vector[N1] x10;
}
parameters{
  real b0;
  real b1;
  real<lower=0> s;
  real<lower=0> sx;
  vector[N] x0;
}  
model{
  for(i in 1:N){
    x[i]~normal(x0[i],sx);
    y[i]~normal(b0+b1*x0[i],s);
  }
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0+b1*x0[i];
    y0[i]=normal_rng(m0[i],s);
  }
  vector[N1] m1;
  vector[N1] x1;
  vector[N1] y1;
  for(i in 1:N1){
    x1[i]=normal_rng(x10[i],sx);
    m1[i]=b0+b1*x10[i];
    y1[i]=normal_rng(m1[i],s);
  }
}
n=20
x0=runif(n,0,20)
x=rnorm(n,x0,2)
y=rnorm(n,x0*2+5,2)
qplot(x,y)

n1=10

#explanatory variable do not has noise
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)

mdl=cmdstan_model('./ex8-0.stan') 
fn(mdl,data)
## Initial log joint probability = -131976 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
## Error evaluating model log probability: Non-finite gradient. 
##       26      -37.5655   0.000105984   0.000189047      0.9582      0.9582       69    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.2 seconds.
##  variable estimate
##    lp__     -37.57
##    b0         5.31
##    b1         1.88
##    s          3.97
##    m0[1]     16.83
##    m0[2]     14.97
##    m0[3]     14.78
##    m0[4]     11.97
##    m0[5]     28.20
##    m0[6]     27.81
##    m0[7]     12.88
##    m0[8]     24.81
##    m0[9]     42.98
##    m0[10]     5.24
##    m0[11]    31.82
##    m0[12]    20.31
##    m0[13]     5.59
##    m0[14]    37.79
##    m0[15]    21.00
##    m0[16]    18.53
##    m0[17]    10.80
##    m0[18]    14.35
##    m0[19]     7.59
##    m0[20]    37.71
##    y0[1]     16.18
##    y0[2]      8.69
##    y0[3]      9.34
##    y0[4]     13.05
##    y0[5]     30.72
##    y0[6]     29.23
##    y0[7]     15.28
##    y0[8]     24.01
##    y0[9]     47.22
##    y0[10]    -1.53
##    y0[11]    23.61
##    y0[12]    25.20
##    y0[13]     4.39
##    y0[14]    32.02
##    y0[15]    17.00
##    y0[16]    25.43
##    y0[17]    13.60
##    y0[18]     6.00
##    y0[19]    10.32
##    y0[20]    36.68
##    m1[1]      5.24
##    m1[2]      9.43
##    m1[3]     13.62
##    m1[4]     17.82
##    m1[5]     22.01
##    m1[6]     26.21
##    m1[7]     30.40
##    m1[8]     34.59
##    m1[9]     38.79
##    m1[10]    42.98
##    y1[1]     10.62
##    y1[2]      8.09
##    y1[3]      9.01
##    y1[4]     19.50
##    y1[5]     17.65
##    y1[6]     25.88
##    y1[7]     29.74
##    y1[8]     40.45
##    y1[9]     45.08
##    y1[10]    42.03
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.4 seconds.
## 
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -37.75 -37.41 1.32 1.03 -40.29 -36.36 1.01      567      945
##    b0       5.14   5.13 1.67 1.65   2.41   7.89 1.01      370      382
##    b1       1.90   1.89 0.17 0.17   1.62   2.18 1.00      509      677
##    s        4.49   4.37 0.83 0.75   3.41   5.93 1.00     1335     1120
##    m0[1]   16.75  16.72 1.06 1.02  15.08  18.45 1.01      566      717
##    m0[2]   14.87  14.85 1.11 1.08  13.08  16.64 1.01      479      624
##    m0[3]   14.69  14.67 1.12 1.09  12.87  16.47 1.01      473      604
##    m0[4]   11.85  11.83 1.25 1.22   9.80  13.85 1.01      409      435
##    m0[5]   28.20  28.19 1.28 1.19  26.10  30.30 1.01     2125     1240
##    m0[6]   27.80  27.78 1.26 1.19  25.73  29.89 1.01     2128     1240
##    m0[7]   12.77  12.75 1.20 1.17  10.79  14.70 1.01      424      460
##    m0[8]   24.79  24.77 1.12 1.06  23.02  26.68 1.01     1865     1214
##    m0[9]   43.09  43.05 2.36 2.28  39.26  46.88 1.00     1054     1219
##    m0[10]   5.07   5.05 1.68 1.66   2.32   7.83 1.01      370      382
##    m0[11]  31.85  31.86 1.50 1.38  29.40  34.22 1.01     1861     1347
##    m0[12]  20.25  20.22 1.02 1.00  18.63  21.98 1.01      932     1054
##    m0[13]   5.43   5.41 1.65 1.63   2.72   8.12 1.01      371      385
##    m0[14]  37.86  37.87 1.94 1.83  34.63  40.95 1.00     1351     1325
##    m0[15]  20.95  20.91 1.03 1.01  19.32  22.68 1.01     1055     1094
##    m0[16]  18.46  18.43 1.03 1.02  16.82  20.15 1.01      700      888
##    m0[17]  10.67  10.66 1.31 1.28   8.51  12.76 1.01      395      390
##    m0[18]  14.25  14.22 1.14 1.10  12.39  16.06 1.01      459      562
##    m0[19]   7.44   7.41 1.51 1.51   4.97   9.87 1.01      375      386
##    m0[20]  37.78  37.78 1.93 1.83  34.56  40.86 1.00     1357     1363
##    y0[1]   16.91  16.85 4.82 4.52   8.99  24.79 1.00     1907     1631
##    y0[2]   14.73  14.88 4.51 4.45   7.13  21.98 1.00     1480     1686
##    y0[3]   14.60  14.57 4.81 4.55   6.58  22.62 1.00     1715     1790
##    y0[4]   11.71  11.72 4.68 4.61   3.86  19.39 1.00     1626     1698
##    y0[5]   28.27  28.39 4.61 4.30  20.77  35.79 1.00     2045     1977
##    y0[6]   27.81  27.81 4.77 4.50  19.99  35.54 1.00     1878     1861
##    y0[7]   12.72  12.63 4.69 4.46   5.14  20.24 1.00     1673     1784
##    y0[8]   24.89  24.88 4.70 4.42  17.11  32.63 1.00     2030     2013
##    y0[9]   42.96  42.85 5.28 5.12  34.64  51.87 1.01     1987     1832
##    y0[10]   4.96   4.98 4.78 4.81  -2.97  12.65 1.00     1330     1671
##    y0[11]  31.81  31.73 4.80 4.79  24.08  39.53 1.00     2042     1821
##    y0[12]  20.29  20.21 4.66 4.45  12.69  27.82 1.00     2150     1864
##    y0[13]   5.56   5.39 4.82 4.77  -2.34  13.42 1.00     1478     1805
##    y0[14]  37.93  37.81 4.90 4.60  30.13  46.13 1.00     1997     1828
##    y0[15]  20.94  20.99 4.76 4.56  13.25  28.68 1.00     1535     1522
##    y0[16]  18.42  18.37 4.75 4.57  10.74  25.81 1.00     1734     1846
##    y0[17]  10.61  10.49 4.73 4.47   3.28  18.39 1.00     1166     1363
##    y0[18]  14.14  14.10 4.60 4.44   6.58  21.68 1.00     1837     1914
##    y0[19]   7.30   7.40 4.73 4.49  -0.37  14.97 1.00     1523     1932
##    y0[20]  37.88  37.86 4.96 4.75  29.86  46.13 1.00     1717     1999
##    m1[1]    5.07   5.05 1.68 1.66   2.32   7.83 1.01      370      382
##    m1[2]    9.29   9.27 1.39 1.39   7.00  11.52 1.01      384      387
##    m1[3]   13.52  13.50 1.17 1.14  11.60  15.40 1.01      440      501
##    m1[4]   17.74  17.70 1.04 1.02  16.11  19.42 1.01      636      790
##    m1[5]   21.97  21.93 1.04 1.02  20.33  23.73 1.01     1291     1147
##    m1[6]   26.19  26.16 1.18 1.10  24.31  28.15 1.00     2051     1213
##    m1[7]   30.42  30.40 1.41 1.30  28.11  32.69 1.01     1996     1328
##    m1[8]   34.64  34.62 1.70 1.59  31.80  37.32 1.01     1617     1328
##    m1[9]   38.87  38.86 2.02 1.91  35.50  42.07 1.00     1265     1363
##    m1[10]  43.09  43.05 2.36 2.28  39.26  46.88 1.00     1054     1219
##    y1[1]    4.97   4.87 4.89 4.67  -2.75  13.06 1.00     1171     1582
##    y1[2]    9.30   9.38 4.76 4.61   1.35  17.06 1.00     1517     1969
##    y1[3]   13.54  13.49 4.79 4.46   5.83  21.18 1.00     1407     1548
##    y1[4]   17.64  17.66 4.64 4.17   9.98  25.38 1.00     1709     1547
##    y1[5]   22.12  22.06 4.63 4.49  14.41  29.82 1.00     1798     1891
##    y1[6]   26.15  26.17 4.79 4.64  18.23  34.27 1.00     1946     1930
##    y1[7]   30.46  30.49 4.98 4.68  22.06  38.43 1.00     1942     1893
##    y1[8]   34.75  34.65 4.79 4.74  27.32  42.63 1.00     2056     2057
##    y1[9]   38.78  38.70 5.12 4.69  30.08  47.08 1.00     1774     1528
##    y1[10]  43.32  43.36 5.17 4.97  34.90  51.55 1.00     1749     1804

#explanatory variable has noise
x10=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x10=x10)

mdl=cmdstan_model('./ex9.stan')
mcmc=goMCMC(mdl,data,wrm=500,smp=1000)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 Iteration:    1 / 1500 [  0%]  (Warmup) 
## Chain 1 Iteration:  501 / 1500 [ 33%]  (Sampling) 
## Chain 2 Iteration:    1 / 1500 [  0%]  (Warmup) 
## Chain 2 Iteration:  501 / 1500 [ 33%]  (Sampling) 
## Chain 3 Iteration:    1 / 1500 [  0%]  (Warmup) 
## Chain 4 Iteration:    1 / 1500 [  0%]  (Warmup) 
## Chain 3 Iteration:  501 / 1500 [ 33%]  (Sampling) 
## Chain 4 Iteration:  501 / 1500 [ 33%]  (Sampling) 
## Chain 1 Iteration: 1500 / 1500 [100%]  (Sampling) 
## Chain 1 finished in 0.5 seconds.
## Chain 2 Iteration: 1500 / 1500 [100%]  (Sampling) 
## Chain 3 Iteration: 1500 / 1500 [100%]  (Sampling) 
## Chain 4 Iteration: 1500 / 1500 [100%]  (Sampling) 
## Chain 2 finished in 0.5 seconds.
## Chain 3 finished in 0.5 seconds.
## Chain 4 finished in 0.5 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.5 seconds.
## Total execution time: 0.7 seconds.
seeMCMC(mcmc,'[m,x]')
##  variable   mean median    sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -42.85 -44.66 11.25 9.48 -58.31 -21.40 1.05       74      155
##    b0       5.36   5.35  1.82 1.81   2.32   8.29 1.00      715      936
##    b1       1.87   1.87  0.18 0.18   1.57   2.17 1.01      733     1149
##    s        3.26   3.34  1.37 1.44   0.95   5.42 1.03      137      108
##    sx       1.50   1.52  0.83 0.94   0.25   2.86 1.05       66      162
##    x0[1]    7.85   7.65  1.65 2.01   5.66  10.58 1.03      197     1456
##    x0[2]    5.46   5.37  1.06 0.88   3.76   7.24 1.01     1568     1743
##    x0[3]    3.78   3.99  1.42 1.48   1.30   5.64 1.02      275      763
##    x0[4]    3.17   3.29  1.12 0.92   1.23   4.85 1.01     1372     1559
##    x0[5]   12.12  12.11  0.99 0.78  10.52  13.78 1.01     2989     2016
##    x0[6]   13.32  13.08  1.50 1.58  11.38  15.92 1.03      246     1205
##    x0[7]    2.66   2.87  1.53 1.58  -0.03   4.64 1.02      255      805
##    x0[8]   10.71  10.60  1.08 0.87   9.05  12.56 1.01     1466     1883
##    x0[9]   20.19  20.09  1.25 0.97  18.29  22.32 1.01     1654     1534
##    x0[10]   0.77   0.63  1.23 1.15  -0.97   2.95 1.01      516      741
##    x0[11]  14.10  14.08  1.05 0.87  12.42  15.84 1.02     2872     1845
##    x0[12]   6.66   6.84  1.40 1.54   4.26   8.51 1.02      232     1116
##    x0[13]   1.33   1.16  1.37 1.45  -0.56   3.74 1.01      294      779
##    x0[14]  16.99  17.05  1.14 0.90  15.13  18.81 1.01     1736     2146
##    x0[15]   8.25   8.26  1.01 0.85   6.64   9.86 1.01     2594     1872
##    x0[16]   5.06   5.26  1.79 2.18   1.94   7.39 1.02      197      736
##    x0[17]   3.26   3.16  1.08 0.90   1.59   5.09 1.02     1266     1574
##    x0[18]   4.10   4.25  1.17 1.06   2.02   5.79 1.01      541     1097
##    x0[19]   1.50   1.43  1.10 0.91  -0.26   3.38 1.01     1237     1235
##    x0[20]  17.96  17.73  1.32 1.08  16.16  20.33 1.01      538     1386
##    m0[1]   20.06  19.85  3.07 3.86  15.53  24.83 1.03      185     1873
##    m0[2]   15.61  15.64  2.01 1.86  12.30  18.87 1.00     1582     2515
##    m0[3]   12.49  12.55  2.58 2.97   8.50  16.52 1.02      284     1397
##    m0[4]   11.35  11.32  2.05 1.86   8.00  14.67 1.01     2435     2354
##    m0[5]   28.04  27.99  1.88 1.71  25.00  31.21 1.00     4183     2241
##    m0[6]   30.26  30.13  2.74 3.14  25.98  34.49 1.02      257     1732
##    m0[7]   10.41  10.53  2.74 3.20   6.17  14.68 1.02      296     1806
##    m0[8]   25.41  25.46  2.01 1.82  22.10  28.64 1.01     1619     1979
##    m0[9]   43.08  43.13  2.36 2.12  39.13  46.84 1.00     2611     1949
##    m0[10]   6.86   6.97  2.41 2.47   2.81  10.54 1.01      417     1956
##    m0[11]  31.73  31.72  2.03 1.89  28.41  34.99 1.01     3918     2699
##    m0[12]  17.86  17.98  2.59 3.00  13.77  21.81 1.02      243     2358
##    m0[13]   7.91   8.05  2.65 2.96   3.51  11.90 1.02      273     1589
##    m0[14]  37.12  37.08  2.14 1.92  33.71  40.68 1.00     2040     2190
##    m0[15]  20.82  20.82  1.91 1.63  17.73  23.92 1.01     3533     2084
##    m0[16]  14.90  15.10  3.27 4.10   9.90  19.68 1.02      202     1650
##    m0[17]  11.51  11.59  2.03 1.87   8.26  14.71 1.01     1263     2033
##    m0[18]  13.07  13.08  2.14 2.18   9.65  16.51 1.01      530     2033
##    m0[19]   8.23   8.32  2.12 1.97   4.69  11.56 1.00     1819     1870
##    m0[20]  38.93  39.10  2.42 2.35  34.79  42.53 1.01      794     2214
##    y0[1]   20.03  20.72  4.72 4.53  11.50  26.46 1.01      449     2597
##    y0[2]   15.56  15.73  4.09 3.45   8.49  22.14 1.01     3522     2430
##    y0[3]   12.43  11.81  4.44 4.03   6.03  20.52 1.01      894     2754
##    y0[4]   11.24  11.10  4.06 3.44   4.73  18.26 1.00     4138     3241
##    y0[5]   28.06  28.02  4.03 3.46  21.46  34.73 1.01     4322     3611
##    y0[6]   30.31  30.84  4.51 4.20  22.35  36.69 1.01      712     2820
##    y0[7]   10.41   9.83  4.52 4.10   3.87  18.51 1.01      784     2821
##    y0[8]   25.48  25.69  4.04 3.41  18.55  31.85 1.00     3597     3159
##    y0[9]   43.05  43.17  4.33 3.73  35.76  50.01 1.01     3471     2935
##    y0[10]   6.79   7.26  4.27 3.72  -0.79  13.15 1.01     1116     2564
##    y0[11]  31.61  31.63  4.09 3.42  24.88  38.29 1.00     3978     3262
##    y0[12]  17.92  17.42  4.45 4.14  11.54  25.72 1.01      785     2850
##    y0[13]   7.87   8.44  4.40 4.00   0.09  14.19 1.00      942     2023
##    y0[14]  37.11  36.90  4.10 3.39  30.50  44.07 1.00     3652     3179
##    y0[15]  20.90  20.88  4.00 3.52  14.47  27.42 1.01     3876     3506
##    y0[16]  14.89  14.31  4.88 4.97   8.17  23.65 1.01      506     1890
##    y0[17]  11.50  11.77  4.00 3.47   4.51  17.84 1.00     3109     3685
##    y0[18]  13.10  12.77  4.10 3.63   6.86  20.20 1.00     2272     2881
##    y0[19]   8.23   8.40  4.12 3.46   1.42  14.86 1.00     3318     3206
##    y0[20]  38.95  39.40  4.31 3.67  31.16  45.41 1.00     1997     3074
##    m1[1]    5.29   5.28  1.83 1.82   2.24   8.22 1.00      714      917
##    m1[2]    9.46   9.48  1.51 1.52   6.91  11.88 1.00      772      896
##    m1[3]   13.63  13.67  1.26 1.24  11.48  15.67 1.00      917     1200
##    m1[4]   17.81  17.83  1.10 1.06  15.91  19.57 1.00     1306     1868
##    m1[5]   21.98  21.98  1.08 1.04  20.15  23.75 1.00     1943     2134
##    m1[6]   26.15  26.16  1.21 1.16  24.13  28.15 1.00     2015     2181
##    m1[7]   30.33  30.35  1.45 1.37  27.97  32.75 1.00     1638     1864
##    m1[8]   34.50  34.49  1.76 1.67  31.66  37.40 1.00     1328     1564
##    m1[9]   38.67  38.69  2.10 1.97  35.32  42.09 1.00     1151     1557
##    m1[10]  42.85  42.88  2.46 2.33  38.87  46.83 1.00     1043     1596
##    x1[1]   -0.06  -0.02  1.69 1.19  -3.03   2.70 1.02     3852     3188
##    x1[2]    2.21   2.20  1.72 1.17  -0.65   4.99 1.02     3935     2477
##    x1[3]    4.42   4.41  1.68 1.17   1.63   7.12 1.02     4216     2514
##    x1[4]    6.65   6.66  1.70 1.17   3.81   9.56 1.01     3647     2613
##    x1[5]    8.87   8.87  1.70 1.15   6.15  11.72 1.02     3522     2349
##    x1[6]   11.11  11.12  1.66 1.15   8.37  13.87 1.01     3853     2517
##    x1[7]   13.34  13.31  1.69 1.20  10.53  16.15 1.02     4113     3224
##    x1[8]   15.60  15.57  1.69 1.19  12.75  18.48 1.02     3930     2392
##    x1[9]   17.76  17.78  1.80 1.23  14.76  20.65 1.02     4020     1965
##    x1[10]  20.03  20.03  1.70 1.15  17.20  22.81 1.02     3728     2793
##    y1[1]    5.33   5.40  3.98 3.51  -1.29  11.94 1.00     2005     2532
##    y1[2]    9.50   9.50  3.86 3.44   3.27  15.98 1.00     2301     2572
##    y1[3]   13.60  13.60  3.78 3.27   7.51  19.80 1.00     3272     2985
##    y1[4]   17.81  17.83  3.69 3.04  11.66  23.85 1.01     3936     3320
##    y1[5]   22.00  22.00  3.69 3.01  15.94  28.10 1.01     3609     2943
##    y1[6]   26.14  26.14  3.74 3.20  19.87  32.17 1.01     3420     2924
##    y1[7]   30.26  30.26  3.75 3.16  24.17  36.53 1.01     3338     2443
##    y1[8]   34.52  34.59  3.94 3.54  28.06  40.92 1.00     2857     2569
##    y1[9]   38.70  38.72  4.05 3.62  32.07  45.26 1.00     2451     2959
##    y1[10]  42.90  42.88  4.26 3.79  35.92  49.97 1.00     2165     3094

y0=mcmc$draws('y0')
smy0=summary(y0)

grid.arrange(
  qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')

m1=mcmc$draws('m1')
smm1=summary(m1)
x1=mcmc$draws('x1')
smx1=summary(x1)
y1=mcmc$draws('y1')
smy1=summary(y1)

xy=tibble(x=smx1$median,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)

qplot(x,y,col=I('red'))+
  geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')

ex10

outlier
y~cauchy(b0+b1*x,s)

ex10.stan

data{
  int N;
  vector[N] x;
  vector[N] y;
  int N1;
  vector[N1] x1;
}
parameters{
  real b0;
  real b1;
  real<lower=0> s;
}
model{
  y~cauchy(b0+b1*x,s);
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0+b1*x[i];
    y0[i]=cauchy_rng(m0[i],s);
  }
  vector[N1] m1;
  vector[N1] y1;
  for(i in 1:N1){
    m1[i]=b0+b1*x1[i];
    y1[i]=cauchy_rng(m1[i],s);
  }
}
n=20
x=runif(n,0,9)
y=rnorm(n,x*2+5,1)
x[1]=3
y[1]=25
qplot(x,y)

n1=10

x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)

mdl=cmdstan_model('./ex8-0.stan') 
fn(mdl,data)
## Initial log joint probability = -20554.5 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       23      -32.9092   0.000287254   0.000301992           1           1       45    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
##  variable estimate
##    lp__     -32.91
##    b0         6.32
##    b1         1.82
##    s          3.14
##    m0[1]     11.78
##    m0[2]     17.34
##    m0[3]     21.92
##    m0[4]      9.34
##    m0[5]     22.01
##    m0[6]     11.87
##    m0[7]      6.80
##    m0[8]      7.67
##    m0[9]     16.43
##    m0[10]    11.21
##    m0[11]     8.21
##    m0[12]    15.36
##    m0[13]    13.96
##    m0[14]    15.47
##    m0[15]    15.20
##    m0[16]    14.40
##    m0[17]    16.83
##    m0[18]    11.48
##    m0[19]     8.46
##    m0[20]    14.82
##    y0[1]     10.33
##    y0[2]     14.17
##    y0[3]     28.26
##    y0[4]     11.59
##    y0[5]     14.15
##    y0[6]     10.50
##    y0[7]      3.89
##    y0[8]      7.90
##    y0[9]     18.40
##    y0[10]    11.94
##    y0[11]     9.42
##    y0[12]    18.08
##    y0[13]    11.17
##    y0[14]    17.45
##    y0[15]    13.38
##    y0[16]    14.20
##    y0[17]    22.26
##    y0[18]     9.26
##    y0[19]     4.56
##    y0[20]     8.61
##    m1[1]      6.80
##    m1[2]      8.49
##    m1[3]     10.18
##    m1[4]     11.87
##    m1[5]     13.56
##    m1[6]     15.25
##    m1[7]     16.94
##    m1[8]     18.63
##    m1[9]     20.32
##    m1[10]    22.01
##    y1[1]     -2.32
##    y1[2]      7.32
##    y1[3]     12.84
##    y1[4]     14.20
##    y1[5]     15.47
##    y1[6]     10.20
##    y1[7]     11.28
##    y1[8]     23.59
##    y1[9]     15.70
##    y1[10]    20.70
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
## 
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -33.33 -33.02 1.25 1.09 -35.78 -31.94 1.01      515      899
##    b0       6.49   6.48 1.54 1.44   3.89   9.10 1.02      257      299
##    b1       1.78   1.78 0.34 0.33   1.23   2.34 1.01      321      465
##    s        3.56   3.46 0.66 0.60   2.67   4.81 1.00     1889     1428
##    m0[1]   11.85  11.81 0.85 0.84  10.50  13.28 1.01      439      734
##    m0[2]   17.30  17.31 1.06 1.00  15.61  19.10 1.00     1521     1291
##    m0[3]   21.79  21.75 1.74 1.67  19.02  24.68 1.01      660     1045
##    m0[4]    9.46   9.42 1.10 1.05   7.66  11.34 1.01      298      299
##    m0[5]   21.88  21.83 1.75 1.68  19.08  24.79 1.01      655     1067
##    m0[6]   11.94  11.90 0.84 0.83  10.60  13.35 1.01      451      745
##    m0[7]    6.96   6.94 1.46 1.37   4.51   9.47 1.02      259      281
##    m0[8]    7.82   7.80 1.33 1.25   5.65  10.11 1.01      266      288
##    m0[9]   16.40  16.42 0.95 0.90  14.86  17.98 1.00     1871     1330
##    m0[10]  11.29  11.25 0.89 0.89   9.87  12.80 1.01      383      594
##    m0[11]   8.35   8.32 1.25 1.18   6.32  10.48 1.01      274      311
##    m0[12]  15.35  15.38 0.86 0.81  13.96  16.80 1.00     1863     1316
##    m0[13]  13.99  13.97 0.79 0.77  12.73  15.31 1.00     1102     1309
##    m0[14]  15.46  15.49 0.87 0.82  14.07  16.94 1.00     1894     1298
##    m0[15]  15.20  15.21 0.85 0.80  13.83  16.64 1.00     1804     1359
##    m0[16]  14.42  14.41 0.80 0.76  13.12  15.80 1.00     1369     1328
##    m0[17]  16.80  16.83 1.00 0.95  15.21  18.48 1.00     1744     1258
##    m0[18]  11.55  11.51 0.87 0.87  10.17  13.02 1.01      406      717
##    m0[19]   8.59   8.57 1.22 1.15   6.58  10.65 1.01      278      311
##    m0[20]  14.82  14.83 0.82 0.78  13.49  16.23 1.00     1603     1353
##    y0[1]   11.81  11.81 3.77 3.55   5.88  17.92 1.00     1720     1798
##    y0[2]   17.34  17.25 3.84 3.54  10.98  23.62 1.00     2106     1695
##    y0[3]   21.73  21.64 3.93 3.81  15.34  28.29 1.00     1794     1688
##    y0[4]    9.44   9.35 3.88 3.54   3.18  16.03 1.00     1435     1686
##    y0[5]   21.92  21.95 4.02 3.87  15.50  28.58 1.00     1515     1701
##    y0[6]   11.85  12.06 3.80 3.49   5.45  17.73 1.00     1844     1875
##    y0[7]    6.83   6.92 3.86 3.59   0.43  12.88 1.00     1173     1769
##    y0[8]    7.87   7.81 3.83 3.71   1.78  14.22 1.00     1363     1801
##    y0[9]   16.38  16.31 3.69 3.63  10.56  22.56 1.00     1894     1661
##    y0[10]  11.32  11.24 3.84 3.77   5.03  17.68 1.00     1981     1894
##    y0[11]   8.33   8.33 3.74 3.54   2.22  14.35 1.00     1246     1677
##    y0[12]  15.26  15.16 3.74 3.58   9.19  21.51 1.00     1994     1863
##    y0[13]  13.96  13.83 3.85 3.80   7.63  20.37 1.00     1930     1860
##    y0[14]  15.50  15.60 3.70 3.51   9.14  21.33 1.00     1906     1689
##    y0[15]  15.20  15.24 3.78 3.62   8.97  21.33 1.00     1788     1638
##    y0[16]  14.57  14.56 3.82 3.70   8.35  20.76 1.00     2059     1842
##    y0[17]  16.74  16.76 3.76 3.49  10.35  22.84 1.00     1860     1892
##    y0[18]  11.56  11.52 3.68 3.56   5.53  17.73 1.00     1862     1675
##    y0[19]   8.51   8.65 3.84 3.67   1.92  14.59 1.00     1152     1798
##    y0[20]  14.86  14.97 3.82 3.70   8.36  21.05 1.00     1854     1941
##    m1[1]    6.96   6.94 1.46 1.37   4.51   9.47 1.02      259      281
##    m1[2]    8.62   8.60 1.21 1.15   6.62  10.67 1.01      279      311
##    m1[3]   10.28  10.24 1.00 0.95   8.69  11.96 1.01      324      361
##    m1[4]   11.93  11.90 0.84 0.83  10.60  13.35 1.01      451      745
##    m1[5]   13.59  13.57 0.79 0.79  12.34  14.90 1.00      891     1208
##    m1[6]   15.25  15.26 0.85 0.80  13.88  16.69 1.00     1827     1315
##    m1[7]   16.91  16.93 1.01 0.96  15.28  18.61 1.00     1706     1347
##    m1[8]   18.56  18.57 1.23 1.16  16.61  20.67 1.00     1069     1118
##    m1[9]   20.22  20.19 1.48 1.41  17.87  22.75 1.00      791     1065
##    m1[10]  21.88  21.83 1.75 1.68  19.08  24.79 1.01      655     1067
##    y1[1]    7.05   7.03 3.87 3.74   0.56  13.27 1.00     1100     1416
##    y1[2]    8.68   8.68 3.74 3.63   2.36  14.62 1.00     1082     1537
##    y1[3]   10.31  10.27 3.68 3.51   4.51  16.42 1.00     1565     1870
##    y1[4]   11.84  11.87 3.75 3.65   5.63  18.01 1.00     1869     1932
##    y1[5]   13.76  13.62 3.74 3.62   7.75  19.90 1.00     1831     1929
##    y1[6]   15.25  15.21 3.78 3.75   9.19  21.44 1.00     2013     1837
##    y1[7]   16.87  16.99 3.77 3.62  10.72  23.24 1.00     1908     1972
##    y1[8]   18.61  18.58 3.78 3.72  12.27  24.66 1.00     1682     1777
##    y1[9]   20.26  20.28 3.88 3.77  13.88  26.66 1.00     1673     1858
##    y1[10]  21.94  21.99 3.97 3.92  15.30  28.12 1.00     1507     1665

mdl=cmdstan_model('./ex10.stan') 
fn(mdl,data)
## Initial log joint probability = -103.485 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       29      -4.37618   4.66905e-05     0.0015828      0.8264      0.8264       40    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
##  variable estimate
##    lp__      -4.38
##    b0         5.08
##    b1         1.94
##    s          0.33
##    m0[1]     10.91
##    m0[2]     16.84
##    m0[3]     21.73
##    m0[4]      8.31
##    m0[5]     21.82
##    m0[6]     11.00
##    m0[7]      5.59
##    m0[8]      6.53
##    m0[9]     15.86
##    m0[10]    10.30
##    m0[11]     7.10
##    m0[12]    14.72
##    m0[13]    13.24
##    m0[14]    14.84
##    m0[15]    14.55
##    m0[16]    13.71
##    m0[17]    16.30
##    m0[18]    10.58
##    m0[19]     7.36
##    m0[20]    14.15
##    y0[1]     15.03
##    y0[2]     16.86
##    y0[3]     22.16
##    y0[4]      7.65
##    y0[5]     21.56
##    y0[6]     10.75
##    y0[7]      5.55
##    y0[8]      6.97
##    y0[9]     22.29
##    y0[10]     9.87
##    y0[11]     0.88
##    y0[12]    14.46
##    y0[13]    11.04
##    y0[14]    14.75
##    y0[15]    13.07
##    y0[16]    13.89
##    y0[17]    15.22
##    y0[18]    10.77
##    y0[19]     7.76
##    y0[20]    13.49
##    m1[1]      5.59
##    m1[2]      7.40
##    m1[3]      9.20
##    m1[4]     11.00
##    m1[5]     12.81
##    m1[6]     14.61
##    m1[7]     16.41
##    m1[8]     18.22
##    m1[9]     20.02
##    m1[10]    21.82
##    y1[1]      7.74
##    y1[2]      6.66
##    y1[3]      7.86
##    y1[4]     -0.33
##    y1[5]    395.05
##    y1[6]     11.11
##    y1[7]     16.77
##    y1[8]     18.20
##    y1[9]     19.83
##    y1[10]    22.63
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
## 
##  variable  mean median     sd  mad     q5   q95 rhat ess_bulk ess_tail
##    lp__   -7.21  -6.77   1.50 1.12 -10.16 -5.67 1.01      528      565
##    b0      5.05   5.04   0.27 0.27   4.62  5.49 1.01      456      262
##    b1      1.95   1.95   0.05 0.05   1.87  2.04 1.01      580      450
##    s       0.39   0.37   0.13 0.12   0.21  0.64 1.01      357      149
##    m0[1]  10.90  10.90   0.15 0.16  10.65 11.13 1.00      522      403
##    m0[2]  16.85  16.85   0.16 0.14  16.59 17.12 1.00     1909     1332
##    m0[3]  21.77  21.75   0.26 0.22  21.37 22.20 1.00     1289     1201
##    m0[4]   8.29   8.29   0.20 0.20   7.98  8.61 1.01      456      288
##    m0[5]  21.86  21.85   0.26 0.23  21.46 22.30 1.00     1275     1188
##    m0[6]  11.00  11.00   0.15 0.16  10.76 11.23 1.00      530      421
##    m0[7]   5.56   5.56   0.26 0.26   5.16  5.98 1.01      455      264
##    m0[8]   6.50   6.49   0.24 0.23   6.13  6.88 1.01      453      268
##    m0[9]  15.88  15.88   0.15 0.13  15.64 16.12 1.00     1907     1452
##    m0[10] 10.28  10.29   0.16 0.17  10.03 10.54 1.00      492      360
##    m0[11]  7.07   7.07   0.23 0.22   6.73  7.44 1.01      453      281
##    m0[12] 14.73  14.74   0.14 0.12  14.50 14.95 1.00     1751     1517
##    m0[13] 13.24  13.24   0.13 0.13  13.02 13.44 1.00      946     1263
##    m0[14] 14.85  14.86   0.14 0.12  14.62 15.07 1.00     1781     1493
##    m0[15] 14.56  14.56   0.13 0.12  14.33 14.77 1.00     1708     1517
##    m0[16] 13.71  13.72   0.13 0.12  13.49 13.92 1.00     1154     1366
##    m0[17] 16.32  16.32   0.15 0.13  16.06 16.57 1.00     1922     1407
##    m0[18] 10.57  10.58   0.16 0.16  10.32 10.81 1.00      504      370
##    m0[19]  7.34   7.34   0.22 0.22   7.00  7.70 1.01      453      280
##    m0[20] 14.15  14.16   0.13 0.12  13.92 14.36 1.00     1365     1402
##    y0[1]   3.06  10.89 336.09 0.55   8.26 13.44 1.00     2015     1932
##    y0[2]  15.39  16.88  56.90 0.59  14.06 19.07 1.00     2108     2055
##    y0[3]  21.68  21.77   5.71 0.61  19.09 24.38 1.00     2032     2014
##    y0[4]   5.20   8.26 150.26 0.60   5.80 10.90 1.00     1653     1960
##    y0[5]  21.95  21.87   5.55 0.67  19.56 24.53 1.00     1980     1890
##    y0[6]  11.16  11.02   8.13 0.58   8.77 13.68 1.00     2051     1822
##    y0[7]   5.65   5.54   9.82 0.67   3.32  7.82 1.00     1926     2056
##    y0[8]   6.11   6.47  14.94 0.64   4.38  9.03 1.00     1664     1554
##    y0[9]  17.21  15.86  91.80 0.62  13.07 18.24 1.00     1926     1877
##    y0[10] 10.49  10.28  19.41 0.61   8.11 12.73 1.00     1891     1932
##    y0[11]  6.69   7.08   9.45 0.65   3.95  9.54 1.00     1580     1879
##    y0[12] 14.80  14.74   9.15 0.57  12.48 16.92 1.00     1792     1841
##    y0[13] 13.52  13.26   9.30 0.59  11.16 15.30 1.00     1674     1981
##    y0[14] 14.49  14.84  23.37 0.56  12.42 17.22 1.00     2157     1958
##    y0[15] 14.55  14.56  14.41 0.57  12.37 17.06 1.01     1900     1742
##    y0[16] 14.40  13.71  24.40 0.58  11.41 16.28 1.00     1947     1709
##    y0[17] 16.44  16.31   6.84 0.59  13.93 19.19 1.00     1926     1931
##    y0[18] 10.90  10.59  16.72 0.59   8.09 13.02 1.00     1767     1873
##    y0[19]  7.52   7.37  13.54 0.65   4.97 10.17 1.00     1930     1927
##    y0[20] 14.70  14.15  26.56 0.58  11.80 16.63 1.00     2103     1931
##    m1[1]   5.56   5.56   0.26 0.26   5.16  5.98 1.01      455      264
##    m1[2]   7.37   7.37   0.22 0.22   7.04  7.73 1.01      453      280
##    m1[3]   9.18   9.18   0.18 0.19   8.90  9.47 1.01      468      310
##    m1[4]  10.99  11.00   0.15 0.16  10.75 11.23 1.00      529      404
##    m1[5]  12.81  12.81   0.13 0.13  12.58 13.01 1.00      812     1281
##    m1[6]  14.62  14.62   0.13 0.12  14.39 14.83 1.00     1724     1505
##    m1[7]  16.43  16.43   0.15 0.13  16.17 16.68 1.00     1922     1407
##    m1[8]  18.24  18.23   0.18 0.16  17.94 18.55 1.00     1807     1449
##    m1[9]  20.05  20.04   0.22 0.19  19.71 20.43 1.00     1612     1281
##    m1[10] 21.86  21.85   0.26 0.23  21.46 22.30 1.00     1275     1188
##    y1[1]   7.44   5.56  81.93 0.66   3.25  8.12 1.00     1466     1614
##    y1[2]   7.57   7.38  19.49 0.62   4.81  9.67 1.00     1670     1612
##    y1[3]   8.82   9.21  18.35 0.58   6.87 11.80 1.00     1964     2015
##    y1[4]  10.86  10.98  23.35 0.60   8.40 13.35 1.00     1978     1814
##    y1[5]  13.04  12.82   9.98 0.60  10.67 15.54 1.00     1812     1853
##    y1[6]  13.24  14.63  50.56 0.58  12.18 17.24 1.00     2103     1931
##    y1[7]  18.83  16.43  79.66 0.58  13.72 19.16 1.00     2008     1767
##    y1[8]  17.83  18.24  16.77 0.62  15.46 20.63 1.00     2016     1958
##    y1[9]  19.79  20.05  27.14 0.60  17.80 23.10 1.00     1886     1659
##    y1[10] 21.96  21.85   5.35 0.67  19.70 24.02 1.00     1889     1971

ex11

censored  

y i=1-N, y~N(m,s)  
  acutualy        ya i=1-Na
  lower censored  yl i=1-Nl, y<L, P(y<L)=cdf(L-m /s)
  upper censored  yu i=1-Nu, y>U, P(y>U)=ccdf(U-m /s)

cdf(z) cumulative normal density function, P((-Infinity,z],z~N(0,1))
ccdf(z) complementary CDF, P([z,Infinity),z~N(0,1))

P(y | x,m,s)=P(ya i=1-Na)* P(yl i=1-Nl)* P(yu i=1-Nu)

ex11.stan

data{
  int N;
  vector[N] x;
  vector[N] y;
  real L;
  int Nl;
  vector[Nl] xl;
  real U;
  int Nu;
  vector[Nu] xu;
  int N1;
  vector[N1] x1;
}
parameters{
  real b0;
  real b1;
  real<lower=0> s;
}
model{
  y~normal(b0+b1*x,s);
  for(i in 1:Nl)
    target+=normal_lcdf(L | b0+b1*xl[i],s);
  for(i in 1:Nu)
    target+=normal_lccdf(U | b0+b1*xu[i],s);
}
generated quantities{
  vector[N] m0;
  vector[N] y0;
  for(i in 1:N){
    m0[i]=b0+b1*x[i];
    y0[i]=normal_rng(m0[i],s);
  }
  vector[N1] m1;
  vector[N1] y1;
  for(i in 1:N1){
    m1[i]=b0+b1*x1[i];
    y1[i]=normal_rng(m1[i],s);
  }
}
n0=20
x=runif(n0,0,9)
y=rnorm(n0,10+3*x,4)
L=15
y[y<L]=L
nl=sum(y==L)
U=30
y[y>U]=U
nu=sum(y==U)
n=n0-nl-nu
qplot(x,y)

xy0=tibble(x=x,y=y)
xya=filter(xy0, y>L & y<U)
xyl=filter(xy0, y==L)
xyu=filter(xy0, y==U)

n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=xya$x,y=xya$y,N1=n1,x1=x1)
mdl=cmdstan_model('./ex8-0.stan') 

mle=mdl$optimize(data=data)
## Initial log joint probability = -6489.61 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       16      -11.0342   0.000321468   0.000423701           1           1       37    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__     -11.03
##    b0        15.27
##    b1         1.16
##    s          2.41
##    m0[1]     17.23
##    m0[2]     19.64
##    m0[3]     22.27
##    m0[4]     21.52
##    m0[5]     24.14
##    m0[6]     20.56
##    m0[7]     20.88
##    m0[8]     17.58
##    y0[1]     18.43
##    y0[2]     23.74
##    y0[3]     18.30
##    y0[4]     20.60
##    y0[5]     22.56
##    y0[6]     20.35
##    y0[7]     22.56
##    y0[8]     18.21
##    m1[1]     15.96
##    m1[2]     16.97
##    m1[3]     17.98
##    m1[4]     18.99
##    m1[5]     20.00
##    m1[6]     21.01
##    m1[7]     22.02
##    m1[8]     23.02
##    m1[9]     24.03
##    m1[10]    25.04
##    y1[1]     13.74
##    y1[2]     15.97
##    y1[3]     22.45
##    y1[4]     19.45
##    y1[5]     19.88
##    y1[6]     20.83
##    y1[7]     17.74
##    y1[8]     25.53
##    y1[9]     25.37
##    y1[10]    30.10
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,'m')
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -12.17 -11.76 1.66 1.38 -15.41 -10.34 1.01      277      317
##    b0      15.11  14.99 3.89 2.89   9.42  21.60 1.02      178      131
##    b1       1.19   1.21 0.80 0.62  -0.10   2.39 1.02      193      180
##    s        3.67   3.27 1.62 1.16   1.97   6.57 1.00      543      571
##    m0[1]   17.12  17.07 2.67 2.03  13.00  21.54 1.03      184      141
##    m0[2]   19.60  19.59 1.53 1.27  17.20  21.97 1.02      318      384
##    m0[3]   22.31  22.33 1.85 1.55  19.36  25.17 1.00      772      704
##    m0[4]   21.53  21.53 1.55 1.33  19.12  23.90 1.01     1199     1013
##    m0[5]   24.22  24.22 2.86 2.33  19.61  28.56 1.01      342      353
##    m0[6]   20.54  20.51 1.40 1.18  18.39  22.71 1.01      836      776
##    m0[7]   20.87  20.83 1.42 1.20  18.72  23.09 1.01     1082      939
##    m0[8]   17.49  17.43 2.47 1.90  13.72  21.60 1.03      187      144
##    y0[1]   17.09  17.14 4.87 3.91   9.91  24.62 1.01      678      603
##    y0[2]   19.56  19.55 4.42 3.55  12.74  26.71 1.00     1306     1191
##    y0[3]   22.28  22.29 4.43 3.59  15.58  29.48 1.00     1693     1674
##    y0[4]   21.36  21.43 4.32 3.57  14.44  27.97 1.01     1933     1675
##    y0[5]   24.04  24.09 5.00 4.04  16.38  31.58 1.00      846      922
##    y0[6]   20.57  20.60 4.27 3.59  13.78  27.24 1.00     1672     1481
##    y0[7]   20.86  20.90 4.26 3.66  14.10  27.42 1.00     1700     1488
##    y0[8]   17.65  17.59 4.85 3.81   9.83  24.82 1.00      644      582
##    m1[1]   15.83  15.75 3.44 2.59  10.73  21.56 1.02      179      135
##    m1[2]   16.86  16.80 2.82 2.14  12.49  21.55 1.03      182      140
##    m1[3]   17.90  17.84 2.24 1.77  14.49  21.63 1.03      191      167
##    m1[4]   18.93  18.92 1.75 1.45  16.30  21.67 1.02      230      233
##    m1[5]   19.97  19.96 1.44 1.22  17.71  22.26 1.02      430      502
##    m1[6]   21.00  20.97 1.44 1.22  18.79  23.26 1.00     1155      989
##    m1[7]   22.04  22.06 1.73 1.46  19.34  24.70 1.00      937      980
##    m1[8]   23.07  23.07 2.22 1.83  19.42  26.47 1.01      511      448
##    m1[9]   24.11  24.11 2.79 2.28  19.64  28.36 1.01      350      354
##    m1[10]  25.15  25.17 3.41 2.72  19.48  30.30 1.01      293      302
##    y1[1]   15.85  15.72 5.28 4.49   7.53  24.24 1.01      369      451
##    y1[2]   16.88  16.82 4.93 4.14   9.07  24.46 1.01      475      489
##    y1[3]   17.83  17.75 4.78 3.79  11.15  25.30 1.01      718     1061
##    y1[4]   18.82  18.79 4.32 3.60  12.17  25.60 1.00     1064     1207
##    y1[5]   19.98  19.88 4.18 3.36  13.42  26.90 1.01     1319     1296
##    y1[6]   21.12  21.08 4.32 3.60  14.59  28.04 1.00     1831     1595
##    y1[7]   22.15  22.05 4.49 3.76  15.22  29.18 1.00     1717     1133
##    y1[8]   23.23  23.13 4.72 3.77  15.78  30.46 1.00     1669     1406
##    y1[9]   24.33  24.23 4.83 3.86  16.69  32.01 1.00      703      906
##    y1[10]  25.17  25.11 5.30 4.24  16.94  33.27 1.01      523      592

y0=mcmc$draws('y0')
smy0=summary(y0)

grid.arrange(
  qplot(xya$y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(xya$y-smy0$median,xlab='obs.-prd.',main='residual')
density(xya$y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')

m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)

xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)

qplot(x,y,col=I('red'))+
  geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')

data=list(N=n,x=xya$x,y=xya$y,
          L=L,Nl=nl,xl=xyl$x,
          U=U,Nu=nu,xu=xyu$x,
          N1=n1,x1=x1)
mdl=cmdstan_model('./ex11.stan') 

mle=mdl$optimize(data=data)
## Rejecting initial value: 
##   Log probability evaluates to log(0), i.e. negative infinity. 
##   Stan can't start sampling from this initial value. 
## Rejecting initial value: 
##   Log probability evaluates to log(0), i.e. negative infinity. 
##   Stan can't start sampling from this initial value. 
## Rejecting initial value: 
##   Log probability evaluates to log(0), i.e. negative infinity. 
##   Stan can't start sampling from this initial value. 
## Rejecting initial value: 
##   Log probability evaluates to log(0), i.e. negative infinity. 
##   Stan can't start sampling from this initial value. 
## Rejecting initial value: 
##   Log probability evaluates to log(0), i.e. negative infinity. 
##   Stan can't start sampling from this initial value. 
## Rejecting initial value: 
##   Log probability evaluates to log(0), i.e. negative infinity. 
##   Stan can't start sampling from this initial value. 
## Rejecting initial value: 
##   Log probability evaluates to log(0), i.e. negative infinity. 
##   Stan can't start sampling from this initial value. 
## Initial log joint probability = -183.679 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
##       26      -23.6838   0.000505193   1.26264e-06           1           1       34    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__     -23.68
##    b0         5.29
##    b1         3.81
##    s          5.74
##    m0[1]     11.70
##    m0[2]     19.60
##    m0[3]     28.25
##    m0[4]     25.76
##    m0[5]     34.35
##    m0[6]     22.61
##    m0[7]     23.68
##    m0[8]     12.86
##    y0[1]     15.34
##    y0[2]     19.96
##    y0[3]     25.82
##    y0[4]     18.53
##    y0[5]     32.44
##    y0[6]     16.79
##    y0[7]     11.90
##    y0[8]     13.00
##    m1[1]      7.56
##    m1[2]     10.87
##    m1[3]     14.17
##    m1[4]     17.48
##    m1[5]     20.78
##    m1[6]     24.09
##    m1[7]     27.39
##    m1[8]     30.70
##    m1[9]     34.00
##    m1[10]    37.31
##    y1[1]      7.78
##    y1[2]     25.63
##    y1[3]     21.57
##    y1[4]     26.13
##    y1[5]     17.69
##    y1[6]     23.14
##    y1[7]     28.77
##    y1[8]     32.05
##    y1[9]     27.18
##    y1[10]    44.19
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,'m',ch=T)
##  variable   mean median    sd   mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -23.99 -23.62  1.54  1.39 -26.92 -22.18 1.01      301      373
##    b0      -1.29   0.36  8.64  7.57 -17.08   9.47 1.02      189      164
##    b1       5.29   4.87  1.93  1.69   3.03   8.92 1.02      211      223
##    s        9.60   8.62  4.26  3.29   4.92  17.78 1.01      328      454
##    m0[1]    7.62   8.63  5.67  5.08  -2.67  14.83 1.02      190      184
##    m0[2]   18.61  18.87  2.98  2.60  13.65  22.94 1.01      339      440
##    m0[3]   30.63  29.79  4.48  3.51  25.15  39.02 1.00      568      572
##    m0[4]   27.17  26.57  3.60  2.83  22.59  33.80 1.01      837      711
##    m0[5]   39.11  37.57  7.14  5.76  30.52  51.88 1.01      341      457
##    m0[6]   22.80  22.66  2.91  2.47  18.38  27.72 1.01      949      691
##    m0[7]   24.27  23.98  3.06  2.51  19.95  29.63 1.01     1037      738
##    m0[8]    9.24  10.18  5.17  4.54  -0.24  15.86 1.02      191      196
##    y0[1]    7.93   9.13 11.56  9.56 -12.07  24.57 1.00      860      596
##    y0[2]   18.60  19.00 10.77  8.62   1.29  35.31 1.00     1926     1184
##    y0[3]   30.53  29.56 11.70  9.41  14.36  50.77 1.00     1182      716
##    y0[4]   27.33  26.77 11.14  9.06  10.79  46.96 1.00     1602     1122
##    y0[5]   38.98  37.43 12.40 10.41  22.03  61.62 1.00      791      915
##    y0[6]   23.12  22.75 11.19  8.65   5.34  41.00 1.01     1906     1233
##    y0[7]   23.83  23.78 11.18  8.84   7.57  41.72 1.00     1822     1342
##    y0[8]    9.06  10.23 11.72  9.75 -10.85  25.86 1.00      793     1071
##    m1[1]    1.87   3.24  7.56  6.78 -11.75  11.41 1.02      188      165
##    m1[2]    6.47   7.52  6.04  5.35  -4.51  14.14 1.02      189      185
##    m1[3]   11.06  11.85  4.63  4.10   2.75  16.98 1.02      195      220
##    m1[4]   15.66  16.20  3.47  3.00   9.62  20.25 1.01      234      345
##    m1[5]   20.25  20.39  2.86  2.51  15.70  24.69 1.01      489      517
##    m1[6]   24.85  24.47  3.15  2.55  20.48  30.68 1.01     1024      700
##    m1[7]   29.44  28.63  4.16  3.22  24.31  37.19 1.01      642      609
##    m1[8]   34.04  32.97  5.50  4.32  27.41  44.09 1.01      433      478
##    m1[9]   38.63  37.16  6.98  5.58  30.23  51.09 1.01      346      472
##    m1[10]  43.23  41.34  8.54  6.93  33.10  58.59 1.01      304      416
##    y1[1]    1.81   3.49 12.90 10.47 -20.83  19.31 1.01      558      605
##    y1[2]    6.20   7.70 12.11  9.79 -15.05  22.50 1.00      706      724
##    y1[3]   11.27  11.71 10.92  9.00  -7.24  27.80 1.01     1070      830
##    y1[4]   15.22  15.61 11.07  8.58  -3.31  32.12 1.00     1534     1187
##    y1[5]   20.33  20.56 10.65  8.70   4.13  37.74 1.00     1673      933
##    y1[6]   24.64  24.42 11.17  8.93   8.45  42.00 1.00     1860     1549
##    y1[7]   29.18  28.60 11.15  9.44  12.88  48.39 1.00     1590     1062
##    y1[8]   34.79  33.42 11.95  9.52  18.17  55.57 1.00     1316     1102
##    y1[9]   38.22  36.74 12.71  9.75  20.87  60.77 1.00     1145     1247
##    y1[10]  43.22  41.08 13.56 10.57  25.90  67.56 1.01      499      514

y0=mcmc$draws('y0')
smy0=summary(y0)

grid.arrange(
  qplot(xya$y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))

par(mfrow=c(1,2))
hist(xya$y-smy0$median,xlab='obs.-prd.',main='residual')
density(xya$y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')

m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)

xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)

qplot(x,y,col=I('red'))+
  geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
  geom_line(aes(x=x,y=m),xy,col='black')

bimodal distribution
mixed normal distribution
data {
  int<lower=0> N;               // データの数
  real y[N];                    // 観測データ
}
parameters {
  real<lower=0, upper=1> theta; // 分布1と分布2の混合比率
  real mu1;                     // 分布1の平均
  real mu2;                     // 分布2の平均
  real<lower=0> sigma1;         // 分布1の標準偏差
  real<lower=0> sigma2;         // 分布2の標準偏差
}
model {
  for (n in 1:N) {
    target += log_mix(theta,
                      normal_lpdf(y[n] | mu1, sigma1),
                      normal_lpdf(y[n] | mu2, sigma2));
  }
}
y0=rnorm(100,0,1)
y1=rnorm(100,-5,1)
y2=rnorm(100,5,1)
y=sample(c(y0,y1,y2),100)
density(y) |> plot()

EM algorithm

library(mclust)

rst=Mclust(y)
summary(rst)
## ---------------------------------------------------- 
## Gaussian finite mixture model fitted by EM algorithm 
## ---------------------------------------------------- 
## 
## Mclust E (univariate, equal variance) model with 3 components: 
## 
##  log-likelihood   n df  BIC  ICL
##            -238 100  6 -504 -505
## 
## Clustering table:
##  1  2  3 
## 37 24 39
rst$parameters
## $pro
## [1] 0.368 0.241 0.391
## 
## $mean
##        1        2        3 
## -4.76533  0.00268  5.18882 
## 
## $variance
## $variance$modelName
## [1] "E"
## 
## $variance$d
## [1] 1
## 
## $variance$G
## [1] 3
## 
## $variance$sigmasq
## [1] 0.819
## 
## 
## $Vinv
## NULL
plot(rst)